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We consider a set of electrostatic problems relevant for determining the real-space structure and 
the ground-state energy of a two-dimensional electron liquid subject to smooth external potentials. 
Three fundamental geometries are investigated: an elongated metallic island, an antidot, and a 
constriction. In the first two cases complete closed-form analytical solutions are obtained, despite 
the absence of rotational or translational symmetries. These solutions govern the shape and size 
of large quantum dots, and also the size of the depletion regions and the density profiles around 
isolated antidots. For the constriction, an exact asymptotical formula for boundary shape is derived 
and arguments are given in favor of its universality. For the cases where the full analytical solution 
cannot be obtained, an approximate method is proposed as an alternative. Its accuracy is verified 
against numerical simulations in a periodic (checkerboard) geometry. 



I. INTRODUCTION 
A. Formulation of the problem 

Studies of physical phenomena in two-dimensional 
metals and thin films often lead to mixed electrostatic 
boundary- value problems. 1 A prototypical example is the 
task of determining the density profile n(r) of a two- 
dimensional (2D) electron liquid in the proximity of ex- 
ternal charges or voltage sources. Such charges and 
sources are used in practical applications to manipulate 
the electron liquid into desired geometrical shapes, e.g., 
quantum dots, narrow wires, constrictions, etc. They are 
also used to intentionally introduce defects, such as an- 
tidotSj^*^ into an otherwise homogeneous system. In 
addition to artificial sources of external potentials, in real 
materials electrons also experience a random potential 
of ubiquitous charged impurities. It appears therefore 
that methods that can tackle the corresponding electro- 
static problems could be of considerable value both for 
applied and for fundamental research. Unfortunately, the 
boundary- value problems are notoriously difficult to solve 
analytically. The goal of the this paper is to show that 
either the complete analytical solutions or their exact 
asymptotics can be found for a number of nontrivial ba- 
sic geometries. 

Our most unexpected unexpected result concerns the 
electrostatics of a constriction. We show in Sec. HVI that 
the long-range Coulomb interaction can cause a signifi- 
cant narrowing of the region occupied by the 2D electron 
liquid in the constriction and that the boundaries of this 
region are described by an interesting nonanalytic func- 
tion, see Eq. (73J). 

The present article extends and g ener alizes the results 
available in the literature, e.g., Refs. lll2U 3 and 6 7 8ll9l ll(l 
and came as an outgrowth of our recent workii on the 
electrostatics of disordered 2D systems. 

We will consider 2D electron systems that are sepa- 
rated from grounded electrodes or other screening bod- 
ies by distances much larger than the interelectron spac- 



ing. Such conditions are realized in semiconductor het- 
erostructures and in field-effect transistors with thick in- 
sulator layers. In these systems electrons interact via 
the l/r Coulomb law and form a metallic liquid if their 
density is not too low. 

The results of the present work are obtained within 
the approximation that the Thomas-Fermi screening ra- 
dius ttf of the electron liquid metallic is vanishingly 
small. This is the correct leading-order approximation 
if all lengthscales of interest exceed ttf- 

In typical semiconductor realizations of 2D systems, 
ttf is of the order of the interelectron separation a c _ c = 
n -1 / 2 , and so our approach is good for studying varia- 
tions of n(r) on lengthscales larger than a e _ e - Note that 
this does not necessarily prohibit us from describing some 
effects that are due to discreteness of electrons. For ex- 
ample, in Sec. m we will be able to calculate the energy 
spacing of the Coulomb blockade peaks because in the 
leading order it is determined by the classical capaci- 
tance. On the other hand, the equations below cannot 
be used, e.g., to study quantum dots with just a few 
electrons. Also, these equations do not apply too close 
to the edges of the metallic regions where n(r) — > and 
so, formally, a _ — > oo. Still, even in these situations 
the solution of the corresponding electrostatic problem 
should provide a valuable insight. 

In the approximation of vanishingly small ttf, the 
electrostatic potential $(r) in the regions occupied by 
the electron liquid is perfectly flat, 

Metal: ra(r) > 0, e$(r) = Me = const, (1) 

where [i e is the electrochemical potential. The rest of the 
2D plane is occupied by the depletion regions (DR) - 
the areas of exponentially small, effectively zero electron 
density that are classically forbidden for the electrons: 

DR: n(r) = 0, e$(r) > fj, e . (2) 

Our goal is to study the conditions that cause DRs to 
appear and their detailed stricture^ 2 To finalize the for- 
mulation of the electrostatic problem we wish to solve for 



2 



this purpose, we need an expression for $ in terms of n. 
We distinguish two cases. 

Case A. — If the number of electrons N e in the system is 
finite, we use 



$(r) = $ ext (r) + t J 



2 / n ( r ') 



d 2 f 



(3) 



where § ext (v) is the external electrostatic potential and 

k is the dielectric constant of the medium. 

Case B. — If we deal with an infinite system with a 

nonzero average electron concentration, we define $ as 

follows: 



$(r) = ~J d 2 r 



,n{v')-a{v') 



(4) 



Here cr(r) represents a spatially nonuniform background 
of opposite charge. In this formulation the primary pa- 
rameter is the average background density n e — (<r(r)), 
and the electrochemical potential in Eqs. (JTJ and (J2J is 
determined by electroneutrality of the system as a whole, 
so that /i e is a function of n e . There is no loss in gen- 
erality in assuming that the background charge density 
a(r) is confined to the same 2D plane as the electrons. In- 
deed, it is easy to see that an arbitrary three-dimensional 
charge density distribution 03 (r, z) creates essentially the 
same electrostatic potential in the z = plane as the fol- 
lowing effective 2D density: 



<r(r) 



d 2 r' 



dz\z\as(r', z) 
27r[(r-r') 2 + z 2 ] 3 / 2 ' 



(5) 



The only difference between the two potentials is the r- 
independent term e 2 C ~ 1 n e , where Cq 1 , given by 



Co 1 



4tt (a 3 \z\ 
« (0-3) 



(6) 



is the inverse geometric capacitance per unit area be- 
tween the 2D layer and the external sources [(. . .) denotes 
the three-dimensional (3D) spatial average] . As the name 
implies, Cq 1 is usually determined by the fixed dimen- 
sions of the structure and the electrodes, and so it is 
almost independent of n e . In this situation, a finite Cq 1 
causes only a trivial linear shift of the electrochemical po- 
tential n e , while the interesting quantity is the deviation 



5/j. = n e ~ e 2 C 1 i 



(7) 



This is one of the quantities we will be calculating be- 
low. We will assume that $> ext (r) and er(r), whichever is 
appropriate, are smooth and bounded functions. 

The major difficulty in solving the above equations 
stems from the mixed boundary conditions and J2J. 
These conditions do not specify where the boundaries of 
the DRs reside or what the potential inside the DRs is. 
They merely state that the DRs may exist and that in- 
side of them e<&(r) must exceed a certain constant value 



fi e . It seems that in the general case, one can make only 
the following two trivial statements. First, in the formu- 
lation J3J) the metallic regions are located in areas where 
the potential energy eQ ex t(r) is sufficiently low, the rest 
of the 2D plane being a DR. Second, in the formula- 
tion (0J, DRs surround negative local minima of <r(r). [If 
a is nonegative everywhere, then the sought ground state 
is simply n(r) = er(r) and the system is free of DRs]. 

As alluded to above, there has been a sizeable amount 
of work devoted to the electrostatics of disordered elec- 
tron systems. Prominent early investigations include 
those of Efros and Shklovski: 13 on the nonlinear screen- 
ing in 3D doped semiconductors and its extensions to 
2jj > i4 i i5 i i6 j n di sor a! ere d. systems DRs typically have some 
irregular shapes, and so there is no any convenient coor- 
dinate system that can be used to take advantage of the 
known techniques of solving integral equations. Thus, the 
random case seems mathematically intractable. Rigorous 
results that have been obtained for such problems are 
limited to certain scaling lawsi^S^ and some asymp- 
totical limitsiii*i& There is a hope however that in reg- 
ular geometries analytical solutions could be easier to 
derive. This is indeed the case. For example, if the exter- 
nal potential ^ ext (x,y) depends only on one coordinate, 
say, x, then the problem can be reduced to an equa- 
tion for just a few parameters, the positions of the DR 
edges^i Once they are determined, the density n{x) at 
all other x can be found from simple analytical formulas, 
see an example below. Similarly, there is a closed-form 
solution in quadratures for the axially-symmetric case, 
&ext = §ext {/) ! provided there is only one DR«i In other 
words, it is known how to find the density profile n(r) 
of a single round droplet or around a circular depletion 
hole. Let us give a few examples. 



B. Examples of exact solutions 

Example 1. — Our first example is a metallic droplet con- 
fined laterally by a parabolic external potential 
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e$ e xt(r) = -U xx r\ U xx > 0. 



(8) 



The density profile of such a droplet is known to be hemi- 
sphericali 



n(r) = ^^U xx \fa? 



(9) 



with the radius a of the droplet related to its electro- 
chemical potential by fx e = U xx a 2 . 

Example 2. — Another instructive example is an isolated 
DR in the form of a perfect circle. As noted above, such 
a depletion hole can form when cr(r) has an isolated neg- 
ative minimum, <7q < 0. Let this minimum be located 
at the origin, r = 0, and let the radius a of the induced 
depletion hole be small enough so that the expansion 



1 



cr(r) = (Tq + o ". 



(10) 
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can be used. Following the formalism of Ref. Q it is easy 
to find that 



-3a /a x 



(11) 



so that n(r) = at r < a. At r > a, n(r) takes the 
formii 



n(r) 



2 2\ a 
- - - arccos - 
z 3 / r 



(12) 



Example 3. — Our last example is the case of a y- 
independent background density 



cr(r) = a(x) = (T + t;VxxX 2 , 



(13) 



which gives rise to the DR in a form of a stripe of width 
2a flanked by the density distributionil 



i(r) = \^J xx \x\\f. 
a 2 = -4a /a xx 



x 2 — a 2 , 



(14) 
(15) 



on the two sides. Note that all the presented exact solu- 
tions, Eqs. |[?JJ, (|12l) . and l|14l) . agree with a well-known 
result that n(r) has a square-root singularity near the 
edge of the DR. Previously, these solutions and their gen- 
eralizations have been used to study the edges of 2D elec- 
tron liquid^ quantum wireSf&S quantum dots^ and an- 
tidotsiM 

In a sense, all the aforementioned examples are one- 
dimensional because § ext (r) [or cr(r)] depend on a sin- 
gle variable. To the best of our knowledge, there are 
no published exact solutions for truly 2D cases, i.e., for 
the geometries where $ ea;t (r) and <r(r) are smooth func- 
tions of position and have no translational or rotational 
symmetries. The present paper is aimed to fill this gap. 
The geometries we consider are as follows. In Sec. [n] 
we study a droplet in a parabolic but not necessarily 
axially-symmctric confining potential. We show that the 
droplet has the elliptic shape, with Eq. © recovered as 
a special case. In Sec. 11111 we derive a formula for the 
density profile around an elliptic depletion hole. This 
formula bridges the limiting cases of Eqs. (|12fl and (|14|l . 
In Sec. IIVI we treat the nonlinear screening problem for 
the saddle-point, which is presumably the most interest- 
ing basic geometry. We derive the asymptotical formula 
for the width of the constriction, Eq. I|73|) , and give argu- 
ments in favor of its universality. In Sec. we examine 
the periodic ("checkerboard") external potential, which 
attains two goals. First, it enables us to study the inter- 
play of the three fundamental geometries (dot, antidot, 
and the saddle-point) examined in Sees. UlTlIVI Second, 
it serves as a testground for an approximate method of 
solving electrostatic problems suggested in our previous 
publication^ 

At the end of each of the following Sections we briefly 
comment on the relevance of the obtained results for vari- 
ous types of experiments. A detailed comparison with the 
available experimental data is deferred for future work. 



II. ELLIPTIC ISLAND 

In this Section we derive the density profile of a metal- 
lic droplet that resides in the external potential 



e<S> ext {r) = hj xx x 2 + ^U yy y 2 , < U xx < U yy . (16) 



If U xx — U yy we must recover Eq. ©; otherwise, if U xx < 
U yyi we expect the droplet to be stretched out in the re- 
direction, along which the confinement is softer. 

The quickest way to obtain the solution for arbi- 
trary U xx and U yy is to use a classic theorem of Frank 
W. Dysoni^S One of the corollaries of this theorem con- 
cerns the 2D charge density distribution 



n(x,y) = n d \ 1 



b 2 ' 



(17) 



which defines an elliptically shaped droplet with semiaxes 
a and b. The Dyson theorem indicates that such a droplet 
creates the in-plane electrostatic potential of the form 



$d(x,y) = {e/K)Tm d 



bdl 



Via 



1 - 



2 + l 2 ){b 2 + l 2 ) 

y 2 



a 2 + l 2 b 2 + I 2 



(18) 



where A is equal to zero inside the droplet and is equal 
to the largest root of the equation 



y 



a 2 + X 2 6 2 + A ; 



= 1 



(19) 



otherwise. This statement can be proved by expanding 
&d( r ) m ellipsoidal harmonics^Si see Sec. lllll and Adp.IaI 
Performing the integration in Eq. I|18(l for the case A = 
0, we get 



zbriA 



K- 



K - Ex 2 E-(l- k 2 )K y 2 



h 2 



b 2 



(20) 



Here K and E are the complete elliptic integrals of the 
first and the second kind, respectively^ 2 evaluated at 



(21) 



(As explained above, we expect a > b). Equations l|16|l 
and (|20[1 indicate that if we choose a, 6, and rid appro- 
priately, we can satisfy the equilibrium condition Q . In- 
deed, we have e<J> = e$> ex t + e&d = Me m the interior of 
the droplet if the following equations hold: 



yy 

2 



(l-k 2 )[K(k d )-E(k d )} 
E(k d ) - (1 - k 2 )K{k d ) ' 
_2_ K(k d ) - E(k d ) 

U. 



XX 

2 



K(k d )k 2 



JT k 2 

V xx™d 



2-K 2 e^{\-k 2 d ){K-E)K' 



(22) 
(23) 
(24) 
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It is easy to see that these equations have a unique so- 
lution for a, b, and n^. Also, from the fact that the 
integrand in Eq. (jlSjl is nonnegative, we conclude ^(r) 
decreases with r, which ensures that the inequality |J2J| is 
also satisfied. Thus, Eq. (|17() is the desired solution. 

Although the parameters of this solution cannot be 
expressed in terms of elementary functions, one can work 
out the limiting cases, which are as follows. Let us define 



particular experimental geometry included. In the first 
approximation, 23 eAV g ~ e 2 /Cd- The analytical asymp- 
totics of this expression are as follows: 



(eAV g f 



7T 2 e 4 



•\/U xx U,, 



6 K 2 

o4 



A', 



Uyy — U x 



~^m 2 ^, U VV ^U XX 



12 k 2 N e 



U x 



(34) 
(35) 



Zd 



Uyy/U xx > 1. 



(25) 



If Zd = 1 — 8z, where Sz -C 1, we have a nearly circular 
droplet with the following parameters: 



He 
2 K 



1 + -Sz + 0{5z 2 ) 



3 

1/4 1/2 



6z + 0(Sz 2 ) 
0{5z 2 



(26) 
(27) 
(28) 



In the limit 8z — > 0, we indeed recover Eq. ©. In the 
opposite limit, Zd ^> 1, we have a strongly elongated 
droplet with parameters 



b 2 = 



2^e 
Uyy 



-Mi 



1 

2z~d 



C = \nz d , (29) 
(30) 
(31) 



Note that by setting z d to infinity and x to zero, we 
obtain the solution for another geometry of basic inter- 
est: an infinite wire in the parabolic confining potential 
e&ext(y) — U yy y 2 /2. The density distribution in a cross- 
section of such a wire is a semicircle^ 



l (y) = 7^2 U yyVb 2 -y 2 , 



(32) 



see Eqs. (fT7jl and (|3T|) . 

Finally, let us calculate the capacitance of the droplet, 
Cd = e 2 dN e /d[i e , where N e total number of electrons in 
the droplet. Integrating n in Eq. (| 1 T|l over the area, we 
get N e = (2Tr/3)n d ab. Using this result, Eqs. p2 |) -pi )l . 
and some simple algebra, we find 



Thus, for the droplet in a parabolic confinement, the 
separation between the Coulomb blockade peaks should 

— 1/3 

scale as JV e , with a coefficient of proportionality that 
depends on the asymmetry of the confining potential. 
This can be tested in experiments where both the size 
and the shape of the quantum dots can be controlled to 
some degree independently and the extensive statistics 
of the Coulomb blockade spacings can be accumulated, 
see, e.g., Ref. l2~i A cautionary note is that one should 
study large dots, N e ^> 1, where neglected here effects of 
disorder, single-particle level spacing, or subtle features 
of the edge structure 2 ^ (see also Ref. l2al are minimal. 



III. ELLIPTIC HOLE 

In this Section we consider the depletion hole that 
forms in a metallic 2D liquid around an isolated neg- 
ative minimum of the background charge density cr(r). 
This problem is relevant for, e.g., determining the den- 
sity profile of the electron liquid around an antidot. 

We assume that the minimum of er(r) is located at 
r = and that the depletion hole is small enough so that 
the expansion 



1 



1 



a(r) = a a + ^a xx x 2 + -a yy y 2 , 0<cr xx <<7 yy , (36) 

can be used (<To < 0). If a xx = a yy , then DR is a circle 
and Eq. (|12|) must be recovered; otherwise, \ia xx < a yy , 
we expect the DR to be elongated in the x-direction, 
along which <r(r) has the slowest growth. Indeed, below 
we show that the DR is the ellipse 



S: £ + £<l, a>b. 
a z b z 



(37) 



Our solution is based on the following series expansion 
of the three-dimensional electrostatic potential $(r, z): 



C d 



K(k d ) 



(33) 



As expected, the capacitance scales linearly with the lin- 
ear size a of the droplet. 

One application of the derived results is the formula for 
the energy separation eAV g of the Coulomb blockade 23 
peaks that would be observed if the droplet is weakly 
connected to external leads. Here V g has the physical 
meaning of the voltage on a gate that controls the size 
of the droplet, with conversion factor appropriate for the 



$ ( r - Z )=EE <i^(A)£M£f>), (38) 

m=0 {p} 

where A, ^t, and v are the ellipsoidal coordinates 

b 2 

< v 2 < k\ = 1 - < /i 2 < 1 < A 2 (39) 

that are related to the original Cartesian coordinates x, 
y, and z by 
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^ = ^AVV, V 2 = k2 J 2 _ kl) {>?-W-kl){kl- V % z * = -^ {X *-l){l-f){l-v>). (40) 

The numerical coefficients in Eq. I|38|l are to be determined, p are certain real numbers (eigenvalues) that depend 
on m and kh, and £^j(£), Ff n (\) are the ellipsoidal harmonics of the first and the second kindsj2£ respectively. The 
definition and some properties of and are reviewed in App. EI 

Below we show that only the following Ef n harmonics 2 ^ are present in the expansion Ij38|l : 



£f(0 = v^e, (4i) 



(0 = Vl^F - Ofc), C± = - ^1 + 2fc 2 T _ ^ + Aki J . (42) 

(the actual values of pi and r± will not be needed). As for F^, the corresponding formulas are rather cumbersome for 
an arbitrary < kh < 1. Fortunately, in this Section we will need primarily their asymptotics at A 2 — > 1, cf. App. lAl 

bf(A) ~ -^=^ +Dl T/>F^i, Dl = -^, (43) 

V ±m „ (l-C±)- 1 , n /^— r n (1 ~ fcg + - 2C ± k 2 )E(k h ) (1 - C±)(l - fc^if (fc/Q 

7 F3 (A) * -VT = T + D±V ^' D±= 2(l-C ± )C ± (k 2 -C ± )(l-k 2 ) ' ™ 



where K and E are again the complete elliptic inte- 
grals m& 

The series (|38|) is designed to satisfy the Laplace equa- 
tion(d 2 /dz 2 + V 2 )$(r, z) = term by term (cf. Refe.Q 
or l27r . Therefore, our task is to demonstrate that with a 
suitable choice of a^, ^(r, z) also satisfies the boundary 
conditions generated by Eqs. CEJ, ©, CH3), and (|37Jl : 

4-$(r,z = ±0) = T27r(e/K)cr(r), reS, (45) 



$(r,z = 0) = 0, r^S. 



(46) 



In the last line we set fi e to zero, for convenience. 

In order to express these boundary conditions in terms 
of A, /-J., and v we note that the points immediately above 
and below the ellipse S [Eq. lpTT)l] correspond to A 2 — > 
1 + 0. Similarly, points immediately above (below) the 
rest of the 2D plane correspond to /i 2 — > 1 — 0. The 
condition will be satisfied if <!> = $!+ $ 2 , where 3>i 
is analytic in A 2 at A 2 = 1, while $2 has a square-root 
singularity, 



$2 



2neaa 



K h 



v/(A 2 -1)(1-m 2 )(1-^ 2 ). (47) 



Expressed as a function of \x and 1/, a takes the form 



v) = (T + ^<Jy y a 2 (n 2 + v 2 - fc 2 ) 



1 a 

2k 



2 
ft 



2, ,2 



Finally, Eq. (|46|l is equivalent to 
*(A,At= l,v) 



(48) 



(49) 



A quick examination of Eqs. (|47() and i|48|) makes the 
above claim that the series l|38[) involves only the ellip- 
soidal harmonics defined by Eqs. I|41|) - (|44() plausible. As- 
suming that this is true, we conclude that on the ellipse 
S, $ must factorize as follows: 



(50) 



where P is a polynomial of the second degree symmetric 
in its two arguments. Since $ must satisfy Eq. (|49|l . it 
has to be of the form 



(l-^) 3/2 



(1- M 2 ) 3 / 2 (1-^ 2 ) 3 / 2 , A = l, (51) 



where $0 is some constant. Combined with the expres- 
sions for that follow from Eqs. and this 
fixes the expansion coefficients a v m to be 



a.. 



-a. 



14 ^1-^+4^ 



(52) 
(53) 



After some more algebra one finds that the boundary 
condition Q47|l can indeed be satisfied if kh is the solution 
of the equation 

cr xx . 2 Mk 2 - l)E(k h ) + (1 - k 2 )K(k h ) 

a yy ( Kh) (l + k 2 )E(k h )-(l-k 2 )K(k h ) ' l ° } 

while a is set by 

2 _ 2a (1 + k 2 )E(k h ) - (1 - k 2 )K{k h ) 



k 2 (l - k 2 )E(k h ) 



(55) 
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For the potential $0 at the center of the DR we get 
2ne k 2 (l-k 2 f/ 2 _ , 



$ 



3 K(l + kl)E-(l-k 2 h )K yv 



(56) 



while at other points on S, Q(r) takes the form stipulated 
by Eq. (EU): 



2 2\ 3 / 2 

*(r)=*o(l- ? -£ 



(r e 5). (57) 



At large distances from the origin, <I>(r, z) is dominated 
by the m = 1, p = p± term in Eq. I|38|) and behaves 
similar to the potential of an electric dipole, 



$(r, z) 



Ph 



k (r 2 + z 2 ) 3 / 2 ' 



47r e 
15 E 



aoab 



(58) 



except the apparent dipole moment ±p^z has the oppo- 
site direction at observation points above and below the 
z = plane. Our formula for ph resembles a well known 
resulliiSiiS for the apparent dipole moment of a round 
hole in a metallic sheet. 

The sought density distribution n(r) outside of S is 
given by 



K $ 

hm - — r . 

2ire z^o \z\ ' 



or, equivalentely, by 



lim 



2nea ^ (X 2 - 1)(1 - p 2 )(l - v 2 ) 
At large r we can use Eq. (|58|1 to obtain 



M / n . Ph 1 



r » a, 



(59) 



(60) 



(61) 



which elucidates how the approach to the perfect screen- 
ing, n(r) — > cr(r), takes place away from the DR. As for 
n(r) near the DR, it is rather complicated, except in two 
cases: (1) a round depletion hole and (2) a slit infinite in 
the x-direction. In the first case the solution is given by 
Eq. i|12[) (see App.^for details); in the second case, re- 
alized at a xx — 0, the solution is given by Eq. I|14f) upon 
replacements x — > y, a — ^ b. 

It is also possible to work out analytically the lowest- 
order corrections to these two limiting cases. Thus, a 
highly asymmetric cr-minimum, Zh = cr xx /a yy <ti 1, gives 
rise to a strongly elongated DR with parameters 



b 2 = 
$0 = 



[1 + Czh + o{z h ) 

4o-q 
a yy 

87re |a | 3/2 
3 K 



lnzr 



1 - 7}£zh + o(Zh) 



(J 



1/2 

IJIJ 



1 



Czh + o{z h ) 



(62) 
(63) 
(64) 



Equation l|15|) is recovered in the limit Zh — * 0. 



On the other hand, if the a-minimum is nearly axially 
symmetric, §z = 1 — Zh 1, then the corresponding DR 
is nearly circular, with the parameters 



3cr 



b 2 = - 



$n = 



3(To 

8 e |a | 3 / 2 



l + -5z + 0{5z 2 ) 
5 



y/3 k (cr xx a y y) 1 /' 1 



l--6z + 0(5z 2 ) 
5 



0{8z 2 ). 



(65) 
(66) 
(67) 



At 5z = 0, Eqs. and iJBSJl reduce to Eq. (fTTft . 

Let us now discuss the effect of an isolated elliptic DR 
on the total energy of the system. It is easy to see that the 
DR causes a positive correction AE to the energy, which 
can be estimated as follows. The unscreened electric field 
that surrounds the DR is of the order of £ ~ 27r(e/K)<7o 
and is concentrated mainly in a volume AV^ ~ a x b x b. 
Hence, 



AE = — [ d 2 r f dz£ 2 (r,z) - — alab 2 
8fj J K 



(68) 



The exact calculation of AE is more convenient to per- 
form in terms of the potential <I>(r) and the total charge 
density en(r) — eer(r), 



AE 



= \e^jd 2 Mr)(l-^ 



y_ 

b 2 



3/2 



16tt 2 e 2 <j 2 ab 2 
105 ~k E{k h )' 



(69) 



The result is in agreement with the preceding estimate. 
[Recall that 1 < E[kh) < n/2, and so E(kh) does not 
have a strong dependence on a/b.] For fixed a xx and 
a yy , we have a, b oc Itrol 1 ^ 2 ! so that AE esc |cro| 7 ^ 2 . 

The energy correction due to DRs affects a number 
of experimentally measurable properties of 2D electron 
systems. One example is the magnetization of quantum 
dots and quantum wires under the quantum Hall effect 
conditions. Although the DRs with n(r) = do not 
appear in that context, a very similar role is played by 
regions of locally depleted topmost Landau leveli^S. 

Another example concerns the DRs induced by random 
impurities, see Fig. ^ Such DRs can strongly influence 
the energy density of macroscopic 2D systems, especially 
at low electron densities. This effect may be important 11 
for the actively studied phenomenon of the 2D metal- 
insulator transition^ 

To show how the obtained formulas can be applied in 
this context let us discuss one experimental technique, 
which is a particular sensitive probe of DRs. This tech- 
nique, pioneered by EisensteinjSi is the measurement of 
the electric field penetration. In order to do this type 
of experimentj22ii22i2i one prepares a structure where the 
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Note that as all the other results in this paper, Eq. I|71() 
is derived in the approximation that ignores non-Hartree 
terms in the energy (kinetic and exchange-correlation en- 
ergies) . This is not a serious omission at low electron den- 
sities n e where the DRs dominate the field penetration. 
However, at higher densities, the non-Hartree terms must 
be included. It may lead to an interesting nonmonotonic 
behavior of AC t ~j (n e ), discussed theoretically in Refs.lilil 
andHH 



IV. SADDLE-POINT 



FIG. 1: The geometry of the field-penetration experiment. 
The 2D layer of interest is sandwiched between the top and 
the bottom metallic gates. The sample contains ionized im- 
purities (dopants) that are shown as randomly scattered plus 
signs. The case of 5-doping is assumed where all the dopants 
reside in a plane parallel to the 2D layer. The impurities cre- 
ate a random electric field that nucleates the depletion regions 
(shown as holes in the probed layer) . These depletion regions 
enhance the penetrating electric field E v that can be detected 
with the help of the bottom gate. 



2D layer is sandwiched between two electrodes that can 
be considered good metals. Once a voltage difference is 
applied between the 2D layer and the top electrode, some 
electric field leaks through the 2D layer and its spatial av- 
erage (E p ) can be measured by monitoring the amount of 
electric charge that has flown into the bottom layer, see 
Fig. ^ It is immediately obvious that DRs must enhance 
(E p ). Using the results of this Section, we are now able 
to calculate how this enhancement of (E p ) is related to 
the concentration and the linear sizes of the DRs. In this 
calculation we will assume that the interlayer distances 
si and S2 (see Fig.^l are large. 

Equation (|58|l gives one way to calculate (E p ). An al- 
ternative and a more general way is to work with the 
energy correction AE and the corresponding contribu- 
tion AC -1 to the inverse capacitance per unit area of 
the device. [AC -1 is the correction to the inverse ge- 
ometric capacitance Cq 1 = Attsi/k, cf. Eq. ©]. For a 
single elliptic DR AE and AC -1 are related by 



ACT 



A d 2 AE 
e= dN 2 



1 d 2 AE Ait 2 



e 2 A da 2 



3nAE(k h )' 



(70) 



where A is the sample area. If the DRs in the 2D layer 
are well separated, their contributions to AC -1 are, in 
the first approximation, additive. The total correction 
to the inverse capacitance, AC^J , characterizes the non- 
ideal nature of our capacitor and therefore the amount of 
field penetration. Indeed, a simple derivation along the 
lines of Ref. |29 leads to the relation 



- — {h p ) - —AC tot - — 

dn e S2 3 K S2 

where TVdr is the DR concentration. 



ab 2 



E{k h ) 



(71) 



In this Section we examine the structure of the 2D elec- 
tron liquid near a saddle-point of the external potential, 



3>ext = \UyyV 2 -\U XX X 2 +0^), U XX ,Uyy>0. (72) 



This is the most interesting fundamental geometry but 
also the most difficult one for the analytical study. 

It is easy to understand that the electron liquid should 
be confined in the y-direction into a strip that has the 
smallest width at x = and widens up as \x\ increases. 
In other words, this structure can be thought of as a lo- 
cal constriction in a 2D electron system. Such structures 
have been fabricated and intensively studied in the past. 
Especially much attention has been devoted to quantum 
point contacts, which are constrictions with the bottle- 
neck width comparable to the interelectron separation.— 
Nevertheless, to the best of out knowledge, the analytic 
solution for the electron density profile n(x, y) and the to- 
tal electrostatic potential &(x, y) around the constriction 
has not been presented. (For a recent numerical work, 
see Ref. 13^) . Instead, the unknown functional forms of 
n and $ have been approximated by expressions chosen 
somewhat arbitrarily. For $(r), a parabolic function^2& 
[as in Eq. (|72|l ] or a parabola with a flat insert^! have 
been used. The boundary of the 2D electron channel 
was modelled variously by combination of wedges^ cir- 
cular arcs^S confocal hyperbolas, or more complicated 

41 

curves^* 

Let y = ±b(x) be the equation for the boundaries of 
the electron liquid in the constriction. Although we have 
not succeeded in solving the electrostatics problem for the 
saddle-point completely, below we show that the correct 
asymptotical behavior of b(x) is given by 



b(x) = cxq\x\ exp —Win 



|x|<x , (73) 



where ao, xq are constants determined by the boundary 
conditions far from the saddle-point. This is perhaps the 
most interesting theoretical result that we achieve in this 
work. 

The exponential factor in Eq. (|73|l is a nontrivial effect 
of the long-range Coulomb interaction (see below). It is 
therefore interesting to do a quick estimate of this factor 
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for a typical experimental setup. As explained in Sec. [I] 
pure electrostatics is not adequate on distances shorter 
than the interparticle separation a G _ G . Therefore, in re- 
ality, Eq. (|75|) applies only as long as b(x), \x\ ^> a _ c ~ 
l/y/n(x, 0). On the other hand, xo can be large. For 
example, using xo ~ and x = a c _ c ~ 50 nrn, which 
are reasonable numbers for GaAs-based point contacts, 
we obtain hx(xo/x) ~ 3.0, so that the exponential factor 
in Eq. I|73ll is approximately 0.2. Thus, the long-range 
interaction can cause a significant narrowing of the bot- 
tleneck of the constriction. 

In the purely electrostatic model (a c _ = 0) we are 
allowed to consider the limit x — > 0. Then Eq. i|73|) 
predicts that the tangent a{x) = b{x)/x of the opening 
angle of the constriction becomes vanishingly small, in- 
dependently of U vv /U xx . This behavior stems from the 
long-range nature of the Coulomb interaction between 
the electrons. With short-range interactions (or with- 
out interactions) electrons fill the constriction up to the 
equipotential contour & ext (x, y) = jjL e ; therefore, the con- 
striction is either bounded by two confocal hyperbolas 
or, when 6(0) — » 0, by the two straight lines. In the 
latter case a = (Uxx/Uyy) 1 ^ 2 , which is finite, in con- 
trast to Eq. H73fl . The physical reason for the difference 
is the ability of the system with long-range interactions 
to modify (screen) the external potential even at points 
r = (x, y) that are outside of the area occupied by the 
electrons. The screening flattens out the confining po- 
tential Q ex t, allowing the electron liquid to spread over 
a larger area compared to the noninteracting limit. The 
nontrivial statement embodied by Eq. (|73|l is that the 
screening is more effective along the longitudinal (x) di- 
rection, so that the extra area occupied by the interacting 
electrons is strongly funneled into the constriction center. 

Let us proceed to the derivation of Eq. (|73|l . We distin- 
guish two cases, which are treated in two separate sub- 
sections below. 



A. Small opening-angle constriction 

First we will treat the limit U xx <^ U yy where & e xt is 
a slow function of a;. In this case the constriction can 
be thought of as an adiabatically narrowing quantum 
wire42 To the lowest order in parameter U xx /U yy , the 
y-dependence of n(x, y) at a given fixed x can be found 
by solving the requisite electrostatic problem assuming 
that $ ex t is x-independent. The solution for n(r) is given 
by Eq. (|3*2|) where 6 should now be understood as a yet 
to be found function of x. The electrostatic potential is 
also easy to find, 

*(r) = Me + \u vy (yVv 2 ^ - b 2 cosh" 1 |) . (74) 



Let us define the one-dimensional (ID) electron density 
q(x) along the wire by 

b 

/l K 
dyn(x,y) = - — UyyX 2 a 2 {x). (75) 

-6 

In view of the last equation, the desired b(x) is directly 
related to q(x). To find the equation for q(x) we substi- 
tute Eq. (|32|l into the boundary condition 0} and take 
y = 0. This yields 

J dx'G{x, x')q{x') = ^U xx x 2 -n e - A${x ). (76) 

-x 

The kernel G(x,x') can be expressed in terms of elliptic 
integrals. At \x — x'\ ^> b(x'), it reduces to the Coulomb 
interaction potential G(x,x') ~ e 2 / n\x — x'\. The term 
A$(xo) in Eq. I|76|l represents the potential created by 
electrons located at points |x| > xq. It is determined by 
the behavior of & e xt at large distances and therefore high 
energies. On the other hand, we are interested mainly in 
the structure of the constriction as small x. Physically, 
we may expect that minor changes in /i e modify b(x) 
near the origin considerably but leave A<j>(x ) virtually 
the same. Therefore, the role of A$(x ) is simply to 
renormalize the electrochemical potential by a constant: 
He -> Me + A$(x ). 

The idea of renormalization proves to be very fruitful 
in the present problem. Indeed, since we are interested 
in the behavior at small distances, we have a freedom 
in choosing the cutoff as long as it exceeds the distance 
of interest. Then /z e should be viewed as a function of 
the cutoff and as such, it must satisfy a certain renor- 
malization group (RG) equation. Similarly, there must 
be an RG equation for U xx and, in fact, for U yy . From 
Eq. I|76l) one finds that to the lowest order in the param- 
eter a(l) <C 1, these equations are 

jjUxx = -U yy a 2 (l), j^yy = -±U yy a 2 (l), (77) 

where I = ]xi\xo/x\, x being the running cutoff. The RG 
flow persists as long as^ |x| ^> b(x). Let us consider the 
most interesting case, 6(0) = 0, where \x\ is always larger 
than 6, so that the RG is able to reach its fixed point. 
What is this fixed point? To the order we are working 
with, Uyy = U yy (x ) + [U xx - U xx (x Q )]/2. This implies 
that for Uyy(xo) ^> U xx (xo) the renormalization of U yy 
can be ignored. Therefore, at the fixed point, where the 
right-hand sides of Eq. I|77|) must vanish, a(l = oo) = 
a(x — 0) = 0, in agreement with Eq. 173|) and statements 
above. 

Let us now derive the complete form of a(x). From 
Eqs. and we find that 

i 

-a 2 (l)[lna(l) + 0(1)] + f dl'a 2 (l') - U **( x °\ (78) 

J Uyy 
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Differentiating both sides with respect to I, we get 



o 



di a 



21n< 



(79) 



which complements Eq. (|77|l and has (|73[1 as the solution. 

It is instructive to rewrite Eq. (|73[) in terms of 
the renormalized U xx [this can be done by integrating 
Eq. 



a 2 (l) = 2 



Uxxjl) 

Uyy 



1 + ln 



Uxx{0) 

U xx (l) 



(80) 



Since U xx (l = oo) = [see Eqs. (J77J and O], Eq. P|l 
describes the decrease of a(x) from its noninteracting 
limit value at x = xo to the asymptotical zero at x = 0. 

The obtained RG fixed point must have a finite basin of 
attraction, presumably a J$ 1, for which Eq. (|73|l must be 
valid. The strong evidence that Eq. j^B} is the universal 
asymptotical law independently of the starting a(xo) is 
furnished by the following analysis of the case a(xo) ^> 1, 
which is the far departure from the found fixed point. 



B. Large opening-angle constriction 



We have a situation where the 2D plane is almost com- 
pletely covered by the metallic liquid, except a narrow de- 
pletion strip that gradually fans out of the origin along 
the j/-axis. This problem is adiabatic with respect to 
coordinate y. To the lowest order in 1/ai, where 



ax 



/r, \l/2 



(83) 



we can solve the system of Eqs. ||TJ, (01, (@J, and JHIJ 
pretending that a is y- independent. The solution is an 
obvious modification of Eqs. (|Tlfl and l|57|): 



n(x,y) 
a 2 (y) 



\f xx \x\ \/x 2 - a 2 (y), 
-Aa(0,y)/a xx , 

2 

+ ~cr xx [a 2 (y) - x 2 } 3/2 . 

6 K 



(84) 
(85) 
(86) 



These equations apply whenever they give real n and $; 
otherwise, n — and $ — /i e . At this level of aproxima- 
tion, the boundaries of the constriction are confocal hy- 
perbolas defined by the equation x 2 = a 2 (y), see Eq. (JSSJ. 
At (To = where the constriction just opens up, these 
hyperbolas become straight lines with a = ai = const, 
while $ acquires the cubic dependence on y: 



If the starting, i.e., long-distance value of a is large, it 
is no longer convenient to parametrize the saddle-point 
by U xx and U yy . Indeed, now the external potential <fr e xt 
is almost completely screened, and so it is rather useless 
in setting up the problem. Instead, as in Sec. IIIII we 
describe the effect of the external sources by an effective 
neutralizing charge density a, so that the total charge 
density in the system is n(r) — c(r), see SecQ] Near the 
origin, <j(r) should be of a saddle-point type, 



f 



f 



a = cr 



2 a vvV 



< (Jyy <C a 



(81) 



In fact, working with <r(r) instead of U(r) also brings us 
closer to the practical side of fabrication of such a large- 
angle constriction. It seems feasible that one can make 
this type of a structure by depositing a static surface 
charge of the form (|81|l nearby the 2D plane but doing it 
with voltage sources (thin metallic gates) may be prob- 
lematic. Although in the former method a xx and a yy in 
Eq. (|5T|) would likely be fixed once the structure is made, 
(To could still be varied by an additional distant gate on 
top of the device. 

As in the case of the small-angle constriction, the rela- 
tion between <ro and /x e , i.e., the top-gate voltage in the 
suggested setup, is determined by behavior at large dis- 
tances. Consequently, in a small interval of fi e of interest 
to us we should have 



(T = (C/e 2 )n e 



(82) 



where C is a constant approximately equal to Co, the 
geometric capacitance per unit area [cf. Eq. ©]. 



(87) 



Let us now show that in a more careful treatment, 
Eqs. I|84|) - (|87|l become invalid at exponentially small x 
and y. It is convenient to assume that at large y, <r(r) 
changes from the decreasing to an increasing function of y 
so that the depletion region terminates at some \y\ = yo- 
In this case the system can be thought of as a 2D metallic 
sheet with an elongated bowtie-shaped hole. Similar to 
the case of a round hole, 44 the 3D electrostatic potential 
at \x\, \z\ 3> a(y) can be sought in the form 



$(r, z) - /i e 



e 

K 



dy'p(y')\z\ 



[x 2 + z 2 + {y~ y') 2 ] 3/2 ' 



-2/0 



which is the lowest-order (dipolar) term in the multi- 
pole expansion compatible with Eq. Q. Since the total 
charge density Ara(r) = n(r) — <r(r) is proportional to the 
discontinuity in d^/dz at z = 0, Eq. I|88() entails 



1 

27 



dy'p{y') 



-yo 



[x 2 + (y - y') 2 ] 3/2 ' 



(89) 



To find p(y) we match the near and far-field asymptotics, 
Eqs. (fH|l and (|S9l . To the lowest order in l/«i, we get 



p(y) = ~(ir/16)a xx a 4 (y). 



(90) 



Next, as in the case of a narrow constriction, we must 
account for the renormalization a — > a + Act of the bare 
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parameters. From Eqs. I|89[l and (|90|l we obtain 



0.06 



dy'[a\y>) - a\y)] ^ 
\y - y'\ 3 + ca 3 (y) ' 



c~l. (91) 



Evaluating this integral and substituting Aer(0, y) into 
Eq. (|85|l . we find a divergent log-correction to a = y/a. 
To facilitate the comparison with Eq. I|73|l , the final result 
can be presented in the form 



b{x) = a l \x\ 1 



~J ln 



(92) 



The most likely behavior consistent with both Eqs. (|T3|) 
and (|92() is as follows. A constriction that appears very 
wide, a = ol\ 3> 1, at large distances, renormalizes first 
into an a ~ 1 structure at x — x r . 



a; exp(-ai/12), 



(93) 



and then into an adiabatically narrowing small opening- 
angle constriction at even smaller x. If so, Eq. (|73[1 is a 
universal asymptotic law. 

Concomitantly, we expect that Eq. (|94|l applies only 
at \y\ ^> x c ; otherwise, it is replaced by 



$(0,y) 



e a 



K 1/2 

u XX 



3/2 



\y\ < x c , 



(94) 



in a nominal agreement with the popular models^Si for 
the electrostatic potential near the constriction. 



C. Numerical results 

In order to verify the outlined analytical theory we 
done a series of numerical simulations. In these simu- 
lations the sought n(x,y), the external potential in the 
form <S> ext = -(l/2)U xx x 2 + (l/2)U yy y 2 + Ux 4 , and the 
Coulomb interaction kernel were discretized on a real- 
space square grid. The units used in the simulations were 
e = k = 1, and the length unit was the size of the grid 
cell. In each run, the energy of the system was mini- 
mized numerically for a trial /j e with the help of MAT- 
LAB™ QUADPROG library function. The value of ^ e 
at which the constriction just opens up was then found by 
minimizing the combination — fi e ) 2 + n 2 at the point 
r = (0,0). To avoid a difficult task of reconstructing 
b(x) from n-data defined on the discrete grid, n(x, 0) was 
studied instead. According to Eq. (|32|) . n(x,0) oc b(x), 
so that the scaling form Eq. I|73[) should equally apply to 
n(x,0). Shown in Fig.^a) is the numerically calculated 



n(x, 0) for the U xx = 1, U yy — 4, the grid size 161 x 81, 
and the system size < 1/2, \y\ < 1/4. On the same 
graph we plot the best fit to Eqs. I|73l) and (|32|l . which 
is obtained for— ao = 1.15, xq = 0.82. Considering that 
for the chosen simulation parameters we are not yet deep 
inside the asymptotic regime a <C 1 [in the interval of 
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FIG. 2: (a) Numerical (dots) and analytical (solid line) results 



for n(x,0) at U xx = 1, U, 



U x 



1, u„ 



0.2, U 



yy 

3.125 



4, U = 2. (b) $(0, y) (dots) at 
The slopes for quadratic and 



cubic dependences are shown for comparison. 



x shown in Fig.^a), a calculated according to Eq. i|32|) 
ranges from 0.5 to about 0.2], the quality of the fit is quite 
acceptable. In other words, we think that our computer 
simulations do support the validity of Eq. I|73|l . 

In Fig. Elb), $(0,y) for the case U xx /U yy = 5 is pre- 
sented. In agreement with Eq. (|87|l . it shows a cubic 
y-dependence at large y and is more consistent with the 
quadratic law of Eq. I|94|l at small y . This again reinforces 
our case for universality of Eq. (|73|l . 

Concluding this Section, we would like to make a few 
brief comments on the relevance of the obtained results 
for experiments. First, our predictions for n(r) and ^(r) 
can be directly verified by a number of currently avail- 
able high-resolution imaging techniques, e.g., near-field 
optical microscopy^ conductance interferometry^£ elec- 
trostatic force microscopy^ or local potentiometry—iS 
Second, it is feasible that the differences between the 
found solution and commonly assumed forms of n and 
$ can also be detected in transport through a quantum 
point contact. The signatures of such deviations and the 
question of how they could be amplified by a suitable 
design of the point contact warrant further study. 



V. CHECKERBOARD 

As an application of the obtained results to a more 
complicated geometry, in this Section we examine the 
2D electron liquid situated on the periodic charge-density 
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FIG. 3: The ground-state density distribution computed nu- 
merically on a 40 x 40 square grid for the symmetric checker- 
board b x = by — 1. Darker areas correspond to higher n(r). 
The area shown in each subplot extends beyond the boundary 
of the unit cell by a half grid cell in each direction, (a) low 
density, n e = 0.0873ni. (b) n e — 0.246ni, which is slightly 
above n p . (c) n e = 0.4707ii, about twice larger than n p (d) 
high density, n e — 0.836ni. The solid lines are n — const 
contours for a set of linearly spaced densities. These densi- 
ties are different in each subplot and are chosen to minimize 
uncertainties in the contour positions that arise due to the 
discreteness of the grid. 



background 



Tlx ( 2ttx 27ry\ 
a = n e + — ^cos — h cos —J , b x >b y . (95) 

We refer to this model as the "checkerboard." It is inter- 
esting because it allows one to study the interplay of the 
three basic building blocks (a droplet, an isolated DR, 
and a saddle-point) that exist for a general <x(r). The 
complete analytical solution of the electrostatic prob- 
lem for the checkerboard geometry remains however un- 
known Instead, we will present a numerical solution 
and will discuss how the results of Sees. can be 

used to understand its structure. We will also derive 
some exact asymptotics and finally, at the end of this 
Section, we will discuss a semianalytic ansatz that re- 
produces many properties of the numerical solution with 
a high accuracy, in particular, its energy as a function of 
n e . 

We start with numerical results, which are shown in 
Fig. [3J These plots represent the distribution of n(r) 
within the unit cell < x < b x , < y < b y that are com- 
puted by a numerical program similar to that described 
in Sec. IIVI In this particular simulation b, 
the unit cell is a square. 

As one can see from Fig.^a), at low density, n e < iii, 
the electron liquid is broken into isolated nearly circular 



by so that 



droplets. The droplets surround the maxima of <r(r) that 
are located at the corners of the unit cell. As n e increases 
at fixed ni, the droplets grow. Their boundaries progres- 
sively deviate from the circular form as they become fun- 
neled towards the nearest saddle-points of c(r), which are 
located at the midpoints of the edges of the unit cell. At 
some density n p (percolation point) the droplets merge. 
In the symmetric checkerboard simulated on 40 x 40 grid, 
this occurs at the average density of 



n p = 0.22ni, b x = b y 



(96) 



From experiments with different grid sizes, we concluded 
that the above value should be close to the percolation 
threshold in the continuum limit but no detailed finite- 
size scaling was attempted. 

Figure |3[b) shows the density profile slightly above n p 
where the continuous path through the electron liquid 
already exists. At n p < n e < ni the most noticeable 
change that takes place as n e continues to increase is the 
contraction of the depletion hole at the center of the unit 
cell, see Figs. |3fc) and (d). Finally, at n e > m (not 
shown) the electron liquid becomes free of the DRs and 
its profile faithfully repeats the background, n(r) = <r(r). 

In the asymmetric checkerboard, b x > b y , the evolu- 
tion of the ground state with increasing n e is similar, 
except that the transition to the global percolation takes 
place in two steps. First, at some density n py droplets 
merge into continuous metallic chains that run parallel 
to the y-axis. Subsequently, at n p > n py , the chains 
become interconnected. This behavior is illustrated in 
Fig. 0| where we display the results of our simulations for 
b x /b y — 2 on the 30 x 60 grid. For this grid size, the two 
aforementioned thresholds were found to be 



' I'D 



= 0.17ni, n p — 0.31ni, 



= 2. 



(97) 



Note that in the asymmetric checkerboard the bound- 
aries of the DRs are elongated along the x-direction. In 
particular, the small droplets at low n e and the small 
DRs at high n e are elliptic in shape. 

Our goal in the rest of this Section is to develop ana- 
lytical approaches that are able to reproduce the above 
numerical findings. 



A. Exact analytical asymptotics 

The structure of the ground state can be determined 
asymptotically exactly in the two limits, the low density 
(n e -c rii) and the high density (ni — n e <C ni). 

We start with the low-density case. Let us split the 
total charge background cr(r) into a part with zero mean, 
cr(r) — n e , and a uniform charge density n e . The former 
produces the electrostatic potential 



$i(r) 
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FIG. 4: Similar to Fig. but for b x = 2, b y = 1. (a) low 
density, n e = 0.0549m. (b) n e — 0.269m, which is between 
n vy and n v . (c) n e — 0.384m, which is above n p . (d) high 
density, n e = 0.897m- 



Function $i(r) has the minima at the corners of the unit 
cell and this is the reason why the metallic droplets that 
form at small n e reside there. Each droplet has the elec- 
tric charge Q — en e b x b y . Consider the droplet centered 
at (0, 0) and denote by & ext (r) the total potential felt by 
the electrons in that droplet due to all the others and the 
uniform n e background. To find $ ea;t (r) we can model 
the other droplets as point charges Q arranged in the 
rectangular lattice. In the leading-order approximation 
in the parameter n e /ri\, within the small area covered by 
the droplet, Q e xt is related to <&i as follows: 



^(r) = $ x (r) + ^-M ( k 



(99) 



Here M(z) is the Madelung constant of the rectangular 
lattice with the unit cell of size lxz -1 . M can be easily 
calculated by the Ewald's method. For example, one 
finds that M(l) = -3.900264920001955. To determine 
the size of the droplet we further notice that <& ex t admits 
the expansion analogous to that in Eq. I|16|) . 



Mo 
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(102) 



Substituting these equations into the formulas of Sec.lTTl 
we find the semiaxes a and b of the droplet to be 



a = b. 
b = a 



3 K(k d )-E{k d )b y 



2ir 2 



1/2 



1/3 / \ 1/3 

"l 



(103) 
(104) 



where kd is the solution of Eq. I|22|) for U xx /U yy specified 
by Eq. (|102fl . The depleted area fraction /dr is related 
to a and b as follows: 



/dr — 1 



nab 
b x b y ' 



(105) 



Using the equation of Sec. [H] we can also calculate the 
corrections to the electrochemical potential and the in- 
verse capacitance, <5/i and AC -1 , respectively, in the 
droplet state: 



AC" 



2tt 2 



h h 2 k 2 



1 1/3 



<5/i = 



3 K{k d )-E{k d ) 
3e 2 



1/3 



2 K 



-AC 



Mo- 



(106) 
(107) 



Equations (|103f> - l|107|) are valid for (bxby) -1 <C n e <C ri\. 
At smaller n e one expects deviations due to the discrete- 
ness of electrons in each droplet. At larger n e there are 
other kind of deviations, from the nonclliptic shape of 
the droplets and their strong mutual interaction. 

Let us now switch to the opposite limit of of high den- 
sity, 5n = n\ — n e <C Hi. In this case we deal with 
small depletion holes that surround the negative min- 
ima of cr(r). Such minima are located at the centers of 
the checkerboard cells, e.g., (b x /2, b v /2). Expanding a(r) 
given by Eq. (|95|l around this point and adhering to the 
notations of Eq. (|36|) , we obtain 



(To = n e — ni = —5n, 

2ir 2 ni (T v „ bl 



b 2 
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(108) 
(109) 



Substituting these expressions into the formulas of 
Sec. IIIII we get the semiaxes of the depletion holes a 
and b to be 



7T 



(l + k 2 )E-(l-k 2 )K 
k 2 h E(k h ) 



1/2 



1/2 



,(110) 

(in) 



where kh is the solution of Eq. Q54JI for (J xx /a yy specified 
by Eq. (|109|l . For the DR area fraction we get 



/dr = 7T 



ab 
b x b y ' 



(112) 
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FIG. 5: The electrochemical potential correction S/j, accord- 
ing to the analytical asymptotics [Eqs. (11U7I and (11141 1 (thin 
lines), numerical simulations (dots), and the trial ansatz 
method (thick solid line) for the checkerboard with the unit 
cell of aspect ratio b x /b y = 2. 



while for S/i and AC 1 we find 



AC = 



4 &; 1- kl I Sn 



3tt b\ E{k h ) 
2 e 1 

5 K 



3/2 



(113) 
(114) 



So far we neglected the interaction among the depletion 
holes. In principle, such an interaction, which is a sub- 
leading correction of a dipole-dipole type, can be included 
perturbatively along the lines of Sec. IIVI However, for 
all b x /by studied in our numerical simulations, it was es- 
timated to be a tiny effect at all densities n e at which the 
approximation of DRs by elliptic holes is still adequate. 
Therefore, we will not discuss such a refinement. 

The comparison between the analytical asymptotics 
and the numerical data for 5[i is shown in Fig. [3] for 
the case b x /b y — 2. As one can see, the droplet pic- 
ture [Eq. (|107fl . the left thin line in Fig. [3] remains ac- 
curate up to n e ~ 0.07ni. The isolated depletion hole 
approximation [Eq. I|114|) . the other thin line in Fig. [3] is 
accurate at n e > 0.4ni. We conclude that our analytical 
asymptotics, which are basically the perturbation theory 
results, indeed work at low and at high expected. 



B. Trial ansatz method 

It is also aparent from Fig. [5] that the derived analyt- 
ical formulas fail at intermediate n e . For example, at 
the percolation threshold, n p ss 0.31n,i, the actual value 
of 5/i is about a factor of two off the nearest analytical 
asymptote. Going to higher orders in perturbation the- 
ory to reduce the discrepancy appears to be cumbersome 
and impracticable. It seems that the quantitatively accu- 
rate description of the ground state of the checkerboard 
model at intermediate n e is currently beyond the reach 
of controlled analytical methods. 



There is however an alternative approach, idea of 
which was introduced in Ref. Strictly speaking, this 
approach is uncontrolled yet it is semi-analytical and as 
we will show below, it reproduces the behavior of Sfj, at all 
n e remarkably well, both for the symmetric and for the 
asymmetric checkerboards. In its simplest implementa- 
tion, this method amounts to adopting the following trial 
ansatz (TA) for the ground state density distribution: 



n a (r) = 6(a - ct D rW cr 2 (r) - a 



DR." 



(115) 



Here 9{x) is the step-function and ctdr is a constant that 
must obey the condition 



dx I dyn a (x 1 y) = n ( 



b x b y . 



(116) 



Clearly, n a (r) is entirely fixed by n e and <r(r) with no 
adjustable parameters. Why choosing the trial state in 
this form? Several reasons can be given. First, it is con- 
sistent with the notion that n(r) is determined primarily 
by the behavior of a at points nearby. (After all, the 
interactions do decay with distance). Second, unless a 
function cr(r) possesses multiple widely different length- 
scales, the behavior of a in the vicinity of a given point 
r is dictated predominantly by the value of cr(r) at the 
same point (e.g., small a tend to be located near min- 
ima, large a — near maxima). Consequently, the purely 
local ansatz n(v) = n a [er(r)] seems reasonable. Third, 
Eq. (|115fl preserves the two asymptotic characteristics of 
the exact solution: a square-root singularity at the edges 
of the metallic regions (cf. Sec.UJ) and the perfect screen- 
ing n — > a at large n (cf. Sec. Illlll . Fourth, one can verify 
that Eq. I|115|l is exact for the DR in the from of an in- 
finite slit [Eq. Q14[l ] and is also rather accurate for the 
round depletion hole [Eq. I|12|) ]. 

Perhaps, the only serious deficiency of the proposed 
TA is the omission of the funncling effect of the saddle 
points. Indeed, according to Eq. (|115fl the boundaries of 
the DRs coincide with the <r(r) = odr = const contour, 
whereas we showed in Sec. IIVI that there are logarithmic 
deviations from such a behavior, and these are noticeable 
in Figs.|21and2] At any rate, the ansatz l|115(l is probably 
the simplest form that one can write down, so it makes 
sense to examine how it performs. Having learned its 
strengths and limitations, one will be in a better position 
to apply this kind of methods in situations where the 
brute force numerical simulations are difficult, such as in 
the models of disordered systems^ 

The implementation of the TA method goes as follows. 
First one selects a reasonably dense set of n e and deter- 
mines the corresponding ctdr by solving Eq. (|116f) on the 
computer. In practice, we did it by approximating the in- 
tegrals in Eq. I|llbj) by a sum over the grid points. Then, 
for each ctdr, one evaluates the total energy energy of the 
corresponding trial state (|115fl . Finally, the electrochem- 
ical potential Sfi is computed by a numerical differentia- 
tion of the total energy with respect to n e . The results of 
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FIG. 6: Comparison between the TA method, analytical 
asymptotics, and numerical simulations for the symmetric 
checkerboard, b x = b y . (a) The correction 8fi [cf. Eq. Q] to 
the electrochemical potential as a function of n e . The mean- 
ing of the lines and the dots is the same as in Fig. 2] (b) The 
scatter plot of n vs. a. Solid curves are the predictions of 
the TA for same densities n e as the dots at the top graph but 
skipping every other n e point for clarity. The leftmost curve 
is for n e = m, the rightmost one (which degenerates into a 
single point) — for n e = 0. The symbols nearby each curve 
are from numerical simulations for the corresponding n e . 

such calculations are shown by the thick lines in Figs. 
andEJa). As one can see, the agreement between the TA 
method results for Sfj, and the corresponding numerical 
datapoints is very good. 

To test the TA method further we can directly compare 
the density distribution n a (r) with the numerically de- 
termined ground state n(r). Such a comparison is shown 
in Fig.|SJb) where we present a scatter plot of n vs. a, for 
the case b x = b y . The spread of the symbols (numerical 
data) with respect to the the solid lines indicates that 
our TA is certainly not exact. However, such a spread is 
not dramatic, and so Eq. (|115f) is a viable approximation, 
especially at low and at high n e . 

One more quantity we can do the comparison for is 
the DR area fraction /dr. As one can see from Fig. 
the TA method performs quite well at all n e , while the 
analytical asymptotics [Eqs. (|105fl and (|112fl ] are obeyed 
in their respective validity domains. 

Finally, let us discuss the estimate of the percolation 
threshold that follows from the TA. According to the 




n e /ni 



FIG. 7: Depleted area fraction according to the analytical 
asymptotics of Eqs. 11051 and 1112H (thin lines), numerical 
simulations (dots), and the TA method (thick line) for the 
symmetric checkerboard. 

TA, the boundaries of the DR are defined by the equa- 
tion <j(r) = (7DFt(?T. e ). Therefore, the percolation oc- 
curs at the average density n* that satisfies the relation 
cdr(«P = n*. Under this condition, the DR bound- 
ary passes simultaneously through all the saddle-points 
in the system. For example, within a single rectangu- 
lar unit cell, the DR has the shape of a rhombus with 
vertices at the midpoints of the cell edges. Solving the 
above equation numerically, we found n* ss 0.31ni. This 
number is independent of z c b = b y /b x because checker- 
boards with different z c b can be mapped onto each other 
by rescaling the coordinate axes. Within the TA, such a 
rescaling does not change the topology of the DRs or the 
average electron density. 

Clearly, the TA is unable to resolve the existence of 
two separate thresholds, n py and n p in the asymmetric 
checherboard. Within the TA, the percolation occurs 
simultaneously in the x- and y-directions. Still, n* is re- 
markably close to the upper (global) percolation thresh- 
old determined numerically for Eq. 
This is also the case at larger z c b, e.g., at z c b — 3 where 
we found n p = 0.32ni. On the other hand, at lower z c b 
the discrepancy grows and reaches its largest relative size 
of about 30% at z c b = 1, see Eq. We believe that 

these discrepancies (one threshold instead of two and the 
value of n p ) originate from the two drawbacks of the TA 
method we mentioned earlier. One is its unability to 
handle widely separate lengthscales, which is the case in 
checkerboards with large z c b- The other is its weakness 
in dealing with the saddle-points. The funneling effect 
of the saddle-points allows the electron droplets to reach 
them sooner as n e increases. Therefore, the continuity of 
the electron liquid is established at a lower n e compared 
to that predicted by the TA. 

One may wonder why the TA method is able to predict 
n p with a much higher accuracjiii (~ 10%) in the case 
of a random a(r). One possible explanation is as follows. 
The funneling effect of the saddle-points that is mishan- 
dled by the TA method is especially pronounced in the 
checkerboard geometry because all the saddle-points have 
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the same value of a, so that the percolation contour has 
to pass through all of them simultaneously. Thus, the 
inaccuracy of our TA is maximized precisely at n = n p . 
In contrast, in the case of a random cr(r) the percolation 
contour passes precisely through the center of a saddle- 
point very rarely, and so the TA works very well. 

We conclude that the TA method is an excellent and 
convenient tool for determining the ground-state energy 
and depleted area fraction but it may be less accurate 
when it comes to more subtle parameters of the real-space 
structure, especially if those are heavily dominated by the 
saddle-points or a hierarchy of multiple lengthscales. 



Acknowledgments 

This work is supported by Chris and Warren Hellman 
Fund and by Alfred P. Sloan Foundation. I am grateful to 
Sen Yang for participation in early stages of this project. 
I thank I. A. Larkin and B. I. Shklovskii for valuable 
comments on the manuscript. 



APPENDIX A: ELLIPSOIDAL HARMONICS 



J 



The ellipsoidal harmonics 2 ^ of the first and the second kinds, and respectively, are defined as the two 

linearly independent solutions of the Lame equation (for A) 



d 2 A Idfdk 



[m(m + I)? - (1 + k 2 h )p]A, f(0 = (e 1)(£ 2 



kl). 



(Al) 



For each m, which has to be a natural number, p can take any of 2m + 1 different values that depend on kh- Functions 
£^(A) and F£(A) at A 2 > 1 are related by 



F£(A) = (2m + l)|i?P(A)| 



dl 



W)\E p m (lW 



(A2) 



At large A, F^X) cx 1/A m+1 . As a rule, F^'s are not expressed in terms of elementary functions. For example, for 
E^ given by Eqs. (glJ-lO, Eq. O leads to the following F£ : 
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where E and F are the elliptic integrals: 22 However, for kh = 0, these formulas simplify to 
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(A6) 
(A7) 
(A8) 



As an application of these formulas, one can derive the 
electrostatic potential <&(r, z) around the elliptic DR dis- 
cussed in Sec. IIIII To do so one needs to substitute 
Eqs. (HI and (HI) for F£, Eqs. gTJ and g2J for Ef nl 
and also Eqs. (|52|l and (|53|l for into the series expan- 
sion (|38l) . Combining such an expression for $(r, z) with 
Eq. i|(j(Jf) one can then, in principle, deduce the formula 



for the density profile n(r) of the electron liquid outside 
the DR. However, this calculation is not presented here 
because for a generic kh the result is rather unilluminat- 
ing. The two notable exceptions are kh = 0, where one 
obtains a circular DR with n(r) given by Eq. 1|12|) . and 
kh = 1 where the DR is an infinite depletion strip and 
Eq. O holds. 
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